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Discrete numerical simulations are performed to study the evolution of the micro-structure and the 
response of a granular packing during successive loading-unloading cycles, consisting of quasi-static 
rotations in the gravity field between opposite inclination angles. We show that internal variables, 
e.g., stress and fabric of the pile, exhibit hysteresis during these cycles due to the exploration 
of different metastable configurations. Interestingly, the hysteretic behaviour of the pile strongly 
depends on the maximal inclination of the cycles, giving evidence of the irreversible modifications 
of the pile state occurring close to the unjamming transition. More specifically, we show that for 
cycles with maximal inclination larger than the repose angle, the weak contact network carries the 
memory of the unjamming transition. These results demonstrate the relevance of a two-phases 
description -strong and weak contact networks- for a granular system, as soon as it has approached 
the unjamming transition. 



I. INTRODUCTION 

The discrete nature of granular materials makes their 
evolutions very complex. In response to an external 
driving force, the macroscopic behaviour of a granular 
system is determined by the intimate interplay between 
the evolving disordered micro-structure, i.e., the arrange- 
ment of the grains and the geometry of their contact net- 
work, and the frictional interactions at the contacts, that 
allow for non trivial force transmissions Q, |^ . Even 
though significant advances have been made in the last 
decade 0:IEIlll3|j the underlying question of a physically- 
based identification of the relevant internal variables and 
of their evolution laws remains open 

A particular issue of recent works on granular media 
concerns the transition from/to a rigid solid- like state 
to/from a flowing fluid- like state, due to its broad inter- 
est in geological and industrial processes and in partic- 
ular to its challenge for condensed matter physics. This 
transition , e.g., the signature of the development of a 
yield stress or flow threshold in a disordered granular 
system, or of an infinite relaxation time compared to 
the actual experimental time-scale, has recently gener- 
ated a flurry of activity [Ulil El IH. Even though 
dense granular materials are athermal, it has recently 
been proposed that the transition could be seen as a 
"u njam ming / j amming transition" at zero temperature 
[Tllll4 [TSi l , motivated by remarkable analogies between 
granular and glassy systems at both macroscopic and mi- 
croscopic scales [la, [13, nil • 

In this picture, the granular system is driven out of 
equilibrium under quasi-static external driving force, in- 
troducing fluctuations, as the system explores different 
packing configurations or metastable states, that could 
be related to an "effective temperature" jlS 0|. Im- 
portant non equilibrium effects of structurally disordered 
granular systems, beside the slow relaxation to equi- 
librium fia l23j |. are : bistability and hysteresis 
[23.I25II26I I as the result of a memory effect due to the in- 



ternal history dependence of systems out of equilibrium 
and the lack of an unique metastable state [J ; and a 
jerky response to the external driving force with impulse- 
like events (local rearrangements of grains or avalanches), 
as the result of the exploration of local minimal energy 
configurations [23, • A related important question is 
whether a diverging leng th-scale exists on the jammed 
side of the transition [29|. 

The well known observation of different characteristic 
angles of the stability of granular systems results from 
such memory effects. When dense granular piles are in- 
clined in the gravity field, there exists a maximal angle of 
stability, or avalanche angle 9a at which the pile starts in- 
evitably to flow, and an angle of repose Or < Oa, defined 
as the slope angle at which the system comes back to 
rest. In the range [6'r,6'a], the system exhibits a bistable 
behaviour where it can be in a jammed -rigid- or an 
unjammed -flowing- state. In that region, the jammed 
state of the granular pile is conditionally stable, e.g., an 
avalanche can be triggered by perturbations of finite am- 
plitude [30I Isij , evidencing metastability of the pile for 
slopes in [9^, 9a\. Under rotation in the gravity field, the 
pile undergoes therefore a subcritical transition. 

In [13, 113, US J the unjamming transition of a 2D cohe- 
sionless granular bed slowly tilted towards the avalanche 
angle 9a was numerically investigated. This transition 
is characterized by a jerky response to the smooth ro- 
tation of the principal stress direction with respect to 
the packing, as attested by the occurrence of local rear- 
rangements of grains. An analysis of the local stresses 
reveals the existence of a significant population of over- 
loaded grains, carrying a shear stress ratio larger than the 
critical threshold of the packing at the unjamming transi- 
tion. This allows to define a coarse-graining length scale, 
or "correlation length" , which increases with the rotation 
up to an angle identified to 9r, where it jumps to a length- 
scale comparable with the thickness of the granular bed. 
This jump can be mapped onto a percolation transition 
of the overloaded grains, and coincides with the onset of a 
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packing dilation, indicating coherent shearing -a macro- 
scopic shear strain across the sample- prior to the unjam- 
ming transition. Further insights in the domain [^^^'a] 
were achieved when analysing the highly frictional con- 
tacts, e.g., the critical contacts where the friction is fully 
mobilized. Indeed the critical contacts are at incipient 
slip and likely give rise to local rearrangements of grains. 
A "correlation length" , based on the multi-scale anal- 
ysis of the critical contacts, exhibits a power-law type 
divergence with an onset that coincides with Or. These 
highly frictional contacts tend to be associated with the 
weak contact network, as defined in 0, [s^ , suggesting a 
two-phases system and a second-order phase transition. 
The emergence of long-range "correlation lengths" within 
the range [6*^,6*0], both in the stress and the frictional 
states, is consistent with experimental measurements of 
both the amplitude of the local perturbation required to 
trigger an avalanche and the increasing scale of the re- 
sponse to a local perturbation, close to the unjamming 
transition [2^ . This suggests that a peculiar regime 
exists before the unjamming transition, namely in the 
range [9r,9a\, where the accumulation of frictional forces 
and resulting local rearrangements of grains leads to co- 
herent shearing and enhances the structural disorder in 
the packing js^ . 

Such modifications of the micro-structure of gran- 
ular systems within the range [0r,^a] at the unjam- 
ming transition should have a signature on non equilib- 
rium effects, such as the system hysteresis during cyclic 
solicitations. In contrast to cyclic shear solicitations 
[tR 113, 13^1 lifil l , we focus here -in the continu- 

ation of [23, ISdl" quasi-static cyclic rotations of 
dense cohesionless granular systems under gravity. De- 
pending on the amplitude of the cyclic solicitations, the 
granular system can approach the unjamming transition, 
without overpass it, so as to keep a solid-like state and a 
quasi-static evolution. The internal history dependence 
or hysteresis can therefore be investigated depending on 
how far from equilibrium the system has been driven, i.e., 
on the distance from the unjamming transition. We re- 
port here results of discrete numerical simulations, based 
on the contact dynamics method, of quasi-static cyclic 
solicitations applied to a 2D cohesionless granular bed. 
These results show that the hysteresis does depend on 
the amplitude of the cycles, and exhibit specific mem- 
ory effects when the system approaches the unjamming 
transition, as a result of the structural evolution. In par- 
ticular, the hysteresis is analysed in terms of both strong 
and weak contact network contributions to the global re- 
sponse. It is found that the weak contact network carries 
most of the signature of the memory effects, confirming 
the relevance of a two-phases description for the unjam- 
ming transition, and extending it to the quasi-static rhe- 
ology, as soon as the granular system has approached the 
transition. 

The paper is organized as follows. The numerical pro- 
cedure and details of the numerical simulations are pre- 
sented in section |n] The section IIIII focuses on critical 



contacts : after discussing previous results during mono- 
tonic external loading, we investigate the evolution of 
the critical contacts during relaxation and cyclic solicita- 
tions. In section llVl the hysteresis of the granular system 
is analysed with respect to the amplitude of the cycles. 
Then the respective contributions of the weak and strong 
contact networks to the hysteresis of the packing are in- 
vestigated in section^ The paper ends with a discussion 
of the results in section IVTI 



II. NUMERICAL PROCEDURES 

Discrete simulations of granular media take into ac- 
count the individual existence of each grain constituting 
the system. The behaviour of the collection of grains is 
entirely driven by the usual equations of motion and the 
contact laws ruling the interactions between the grains. 
The contact dynamics method, applied for this work, 
deals directly with infinitely stiff contact laws, assum- 
ing that grains interact through hard-core repulsion and 
non-smooth Coulomb friction only • This implies that 
the contact force between two grains is non-zero only if 
the latter are touching. Once a contact is created be- 
tween two grains, the latter cannot get closer, so that 
any normal relative motion is only repulsive : the grains 
are perfectly rigid and cohesionless. The Coulomb fric- 
tional law consists of an inequality between the tangen- 
tial and normal forces at contact, referred to as T and 
N respectively : ^(xN < T < ifiN, where /i is the mi- 
croscopic coefficient of friction, and sets a threshold for 
a contact to slip. Accordingly, this controls the frictional 
dissipation at the contact scale in the case T = zLfxN 
when tangential slip motion is possible according to the 
immediate environment of the contact, by contrast with 
the case ^ij,N < T < zL^N and no slip motion is allowed. 
In the present work, we are interested in the quasi-static 
evolution of the packing and concentrate on micro-plastic 
rearrangements only, therefore we consider purely dissi- 
pative collisions, and set the coefficient of restitution to 
zero. 

The performed simulations are bi-dimensional. We 
consider packings of 4000 circular grains, randomly de- 
posited in a rectangular box. The grains diameters are 
uniformly distributed in the interval [dmm', dmax] with 
dmax/dmin — 1-5. The polydispcrsity, measured as 
(dmax — drain) /i2d) — 20% with d the mean disk diame- 
ter, prevents long-range crystal-like ordering in the pack- 
ing by enhancing structural disorder. Each granular bed 
shows a nearly flat surface parallel to the bottom of the 
box, with a slope 6 with respect to the horizontal direc- 
tion ; its thickness is ~ SSd, and its length is ~ 120d. The 
granular packings have an initial solid fraction C ~ 0.785 
and an initial coordination number Z ~ 3.46. The micro- 
scopic coefficient of friction (/i = 0.5) is the same between 
all the grains and between the grains and the walls of the 
box. Figure ^ shows a typical granular pile inclined in 
the gravity field g at an angle 6 = 17°, with the xy- frame 
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FIG. 1: A typical simulated granular pile of 4000 disks whose 
diameters are polydispersed at 20%, inclined at an angle 9 — 
17° in the gravity field g with the xy-frame linked to the 
bottom of the box and the unity vector normal to the free 
surface ns- 




is resumed. Two successive cycles are performed. This 
solicitation is applied on 50 granular piles differing in 
the initial disordered micro-structure, i.e., the grains ar- 
rangement and the contact network geometry, and show- 
ing characteristic angles 9r — 15° and 9a — 21°. In the 
following, averages are taken over all 50 simulations. The 
granular beds and the data describing the first loading 
stage, are the actual data of Staron et al. ^32] . 

We study the influence of the cycle amplitude and of 
the distance from the unjamming transition on the evo- 
lution of the pile state by considering two values of the 
maximal inclination angle O^ev , namely by applying small 
and large cycles. The smallest value of Orev is chosen so 
that the small cycle keeps the pile out of the coherent 
shear regime identified prior to the unjamming transi- 
tion : 6*^611 — 10° < 9r. On the contrary, the largest 
value of Orev is chosen in the range of slopes [9r',9a]^ 
but below the minimal avalanche angle observed among 
the 50 realizations to prevent the pile from avalanching : 
Orev — 18° e [Or] Oa\. These two kinds of cycles are repre- 
sented in FigureEl where the slope of the pile is plotted 
against the time t for both small (in gray) and large (in 
black) cycles. Continuous lines represent loading stages, 
and dashed lines represent unloading ones. The horizon- 
tal lines indicate the repose and avalanche angles Or and 
Oa, with a thickness equal to the size of the fluctuations. 
A diamond indicates the initial state of the granular pile 
just before the cycles start, that will be reported on all 
the next figures. 



III. CRITICAL CONTACTS CLUSTERING 



FIG. 2: Temporal plot of the surface slope 6(t) of the granular 
pile during the loading-unloading solicitations for Orev = 10° 
(in gray) and Orev = 18° (in black). The horizontal lines, 
with a thickness equal to the size of the fluctuations among 
50 realizations, correspond to the repose and avalanche angles 
Or and Oa. The initial state of the pile just before the cycles 
start is represented by a diamond. 



linked to the bottom of the box and the unity vector nor- 
mal to the free surface ns- In practice, instead of rotating 
the granular pile in the constant gravity field, the direc- 
tion of g with respect to the pile is rotated to simulate 
tilting. 

The cyclic solicitations consist of rotations at a con- 
stant rate (±0.001° per time step), the rotation being 
positive when clockwise. A positive rotation is first ap- 
plied to newly generated granular piles so that the slope 
of the free surface increases from = 0° to a maximal 
inclination angle Orevi corresponding to a loading stage. 
After the inclination Orev is reached, the direction of rota- 
tion is reversed, and the slope decreases back to 0°, corre- 
sponding to an unloading stage. Rotation is maintained 
so that the slope of the pile evolves from towards the 
opposite maximal inclination angle —Orev The direction 
of rotation is again reversed, and the cyclic solicitation 
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FIG. 3: Density of critical contacts v averaged over 50 simu- 
lations (a) during the first loading starting from S ~ 0° and 
(b) during one large cycle, starting from a pile inclined at 
—Orev = —18°, when varying the slope inclination towards 
the avalanche angle. 

As discussed in the case of slow monotonic external 
loading of a granular slope in [23, , the coherent 
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shear regime identified prior to the unjamming transi- 
tion corresponds to the structuration in the packing of 
highly frictional contacts, e.g., contacts where forces have 
reached the Coulomb frictional threshold, i.e., T = ±fJ,N, 
referred to as critical contacts. An interesting issue would 
be the evolution of the critical contacts during cyclic so- 
licitations. 

We discuss first the robustness of previous results on 
critical contacts ll^,!!!] against the pile preparation. Fig- 
ure 01 displays the evolution of the density of critical con- 
tacts ly, defined as the proportion of critical contacts to 
the total number of contacts in a granular packing, av- 
eraged over 50 piles, when varying the pile inclination 
6 towards 6a, from 0° (a) and from —Orev — —18° {h). 
At the beginning of the solicitation, both evolutions of v 
show a sudden drop to a minimum value indicating the 
instantaneous loss of critical contacts. After this first 
stage and in spite of the different initial conditions, the 
behaviour of the critical contacts when approaching the 
unjamming transition is approximately the same : v in- 
creases up to a maximal value — 27]. 

The loss of critical contacts occurring at the transition 
from loading to unloading in Figure also observed in 
[so*!, indicates that the averaged frictional state of the 
packing related to past solicitations is immediately mod- 
ified. It is related to a systematic grains rearrangement, 
during which critical contacts leave the Coulomb fric- 
tional threshold. Note that the amount of critical con- 
tacts lost at the transition is not larger than the typical 
fluctuations around the averaged critical contacts density 
y during the cycle. In fact, the transition from loading 
to unloading imposed at Orev synchronises the grains re- 
arrangements for the various realizations, hence a larger 
response on average than in the other parts of the cycle. 
The effect of the reversal on a single realization remains 
an open issue, which is beyond the present averaged anal- 
ysis. 

We now observe the relaxation to a static packing in 
terms of critical contacts by suddenly stopping the ro- 
tation. Figure 01 shows the evolution of as a func- 
tion of time after the rotation was stopped at t = 
in 5 granular piles inclined at different angles {9 — 
3°, 11°, 16°, 18°, 20°). The state of the pile before the re- 
laxation starts is represented by a circle. Whatever the 
density of critical contacts and the inclination when the 
packing starts to relax, v vanishes after a complex tran- 
sient dynamics. The time needed for critical contacts 
to relax is larger for slopes close to Oa : this confirms 
the complex relaxation towards equilibrium of a granular 
packing (12 . i22^ .23.] , and this is consistent with observa- 
tions of long- range c orrelations close to the unjamming 
transition [211271 Isotlslls^l. The relaxation of the criti- 
cal contacts for rigorous static conditions shows that the 
critical state is a transient state. This can be understood 
by the fact that the contacts can not sustain the sliding 
condition once the loading is stopped. In this sense, the 
density of critical contacts does not characterize the pile 
plastic state, because is a dynamic response function to 
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FIG. 4: Relaxation of the density of critical contacts v as & 
function of time in 5 granular piles inclined at different angles 
[9 = 3°, 11°, 16°, 18°, 20°) after the rotation was stopped at 
t = 0. The state of the pile just before the relaxation starts 
is represented by a circle. 



the actual loading of the pile. Also other variables could 
be more representative of the plastic state of the pile, 
e.g., the averaged cumulative slip dislocation at critical 
contacts. However critical contacts remain active during 
quasi-static loading, and it was shown that the evolu- 
tion of V for one realization |23 becomes more and more 
intermittent, exhibiting successive frictional demobiliza- 
tion and remobilization periods, when the rotation rate 
decreases. 
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FIG. 5: Density of critical contacts v averaged over 50 simu- 
lations, slide-averaged over l°-intervals as a function of the 
slope angle 6 during the successive cyclic solicitations for 

Orev 18 . 



Figure shows the evolution of as a function of 
the pile slope during loading and unloading process for 
6rev — 18°. The density of critical contacts at an an- 
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FIG. 6: Size of the correlated clusters of critical contacts rc, 
slide- averaged over l°-intervals as a function of the slope angle 
9 during the successive cyclic solicitations for 6rev = 18°. 

gle is smaller during unloading than during loading, as 
attested by the hysteretic loops. Figure displays the 
evolution as a function of 6 during the large cycles, of 
the "correlation length" rc, defined as the mean size of 
clusters such that locally i/ = Vc- During each loading 
stage, rc increases sharply above 0r and is expected to 
reach the size of the pile at the unjamming transition, 
as shown in the case of monotonic external loading [2^ . 
At the transition from loading to unloading, a signature 
of the loss of critical contacts is visible : rc immediately 
reduces to a few grain diameters. There is no more evi- 
dence of correlation within the granular pile in terms of 
critical contacts upon reversal of the solicitation. Note 
that after the first loading, i/ as well as rc shows the same 
cyclic evolution, indicating the same behaviour of criti- 
cal contacts during the two successive cycles, especially 
when approaching the unjamming transition. The reduc- 
tion of frictional mobilization after the reversal, induces a 
strong hysteresis in both the evolution of v and during 
loading and unloading cycles. Is this hysteresis apparent 
on other global characteristics of the granular pile? 

The critical contacts play a significant role when ap- 
proaching the unjamming transition, since they likely 
give rise to micro-plastic events. Accordingly, one may 
wonder whether the response of the pile is aed by the hys- 
teretic evolution of critical contacts observed in the range 
of slopes [9r;9a]- Also one would like to know whether 
the possible effect of the observed hysteresis is localized 
in the domain [6r',9a\ or conditions the response of the 
pile along the full cycle of solicitation. In the follow- 
ing, we investigate the hysteretic response of the pile to 
cyclic solicitations and the effect on this response of hav- 
ing approached the unjamming transition. To do so, we 
consider two different amplitudes of loading-unloading 
cycles, namely two values of the reversal angle 0rev, a-nd 
we compare the response of the granular pile in terms 
of stress, strain and evolution of the contact network for 



9rev = 10° and Orcv = 18°. 

IV. HYSTERETIC PHENOMENA DURING 
LOADING AND UNLOADING SOLICITATIONS 

During cyclic tiltings, the induced stresses as well as 
the resulting deformations of the granular pile are inves- 
tigated. We start with the comparison of the global stress 
state during loading and unloading stages. 

A. Stress state of the granular pile 

The maximal value of the shear to confinement ra- 
tio, the so-called inertial number in jl3j is defined as 
Imax = ^maxdy/^, with jrnax, thc maximal shear strain 
rate evaluated during coherent shearing in the range of 
slopes [9r',9a], P, the volumetric mass of the grains and 
P, the typical pressure induced by gravity. In the present 
case, Imax — 10~^ shows that the pile evolution is actu- 
ally quasi-static. As a result, static stresses transmitted 
via contacts are much larger than kinematic ones due to 
momentum transport during collisions, so that the stress 
tensor of the granular pile cr in the xy-ba,me is evaluated 
as follows from |42| : 

1 

'^^^ = vT.ftif, (1) 

a—l 

where i and j denote the space dimensions and a the 
contact. / is the contact force, I is the vector normal to 
the contact surface, and n'^ is the total number of contacts 
in the representative element volume V . The tangential 
and the normal stresses <tt and <7n along the direction 
of the free surface are calculated : 

crjv = II {(Tns)ns \\, (2) 
= II crns - (JNTis II, (3) 

where ns is the unity vector normal to the free surface. 
The direction of the eigenvector of cr corresponding to 
the largest eigenvalue is the principal stress direction, 
whose angle with the direction ns is VP. The principal 
stress direction corresponds to the direction where stress 
is purely compressive. 

Figure |3i displays the shear stress ratio otIoj^ and 
Figure [71) the evolution of the principal stress direction 
* during smaU {drev = 10°) and large {B^ev = 18°) 
loading- unloading cycles. Both ut/un and 4* exhibit 
hysteretic evolutions related to irreversible modifications 
of the force transmissions at contacts, but the hysteresis 
loops are of rather small amplitude and independent of 

rev • 

During the first loading, ut/ctn — atan(0) in agree- 
ment with the predictions of continuous media mechanics 
for an infinite slope inclined in the gravity field at static 
equilibrium, except that a is not exactly equal to one : 
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FIG. 7: Shear stress ratio ctt/o-jv (a) and principal stress di- 
rection ^ (b) as a function of the pile slope 9 during cyclic 
rotations for small {9rev ~ 10°) and large {Orev ~ 18°) ampli- 
tudes. The vertical lines correspond to the value of —6r and 
8,.. The initial state is represented by a diamond. 



a — 0.8. It can be shown that this difference is to be 
attributed to finite size and walls effects (33] ■ During the 
following cycles, ar/cTN — ata,n(0) + e{aT /crN)res, where 
e = — 1 (rcsp. e = +1) during loading (resp. unloading) 
and a — 0.8. A residual stress ratio (ut / <JN)res of the or- 
der of 5% of the typical shear stress ratio is recorded for 
~ 0°. The principal stress direction ^E* tends to ±45° for 
large inclinations, exhibiting a rather smooth evolution. 
For hydrostatic equihbrium, i.e., equal normal stresses : 
CTra = ^yyi it would be expected that ^E" = —45° (resp. 
= +45°) for 6* > (resp. 6* < 0). 

The stress state evolution of the granular pile is plotted 
on the figures [3 for the two successive cycles : the plots 
perfectly collapse after the first loading. In the following, 
the two successive cycles are reported on all the figures. 

Altogether both the shear stress ratio ar/cTN and the 
principal stress direction VE" are controled by the inclina- 
tion 9 (and its temporal derivative, as evidenced by the 
weak hysteresis) , whatever the solicitations path and the 



history of past cycles. 



B. Deformation of the granular pile 

The exploration of different metastable configurations 
results in the apparition of critical contacts in the pack- 
ing, that induces local rearrangements of grains (see sec- 
tion Hill and 113). The effect of these rearrangements at 
the scale of the pile is apparent when plotting the varia- 
tions of the volumetric strain ey = {V — Vo)/Vo, where 
V is the volume of the pile and Vb its initial volume, as 
a function of the slope of the pile during successive small 
and large cycles (Figure IS)). At the first sight, we notice 
that, even if the initial solid fraction C and the stress 
state evolution of the pile are the same for the two cyclic 
amplitudes (Fig. 01 , the volumetric strain evolution does 
depend strongly on Orev 
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FIG. 8: Volumetric strain ev as a function of the pile slope 
during small and large cycles (same symbols as figure|7|l. The 
two lines are indications of the compaction rate dey/dt for 

Orev = 10° and Orev = 18°. 

Irrespective of the value of Orev , cycles produce an over- 
all densification of the pile of the order of 10"'^, corre- 
sponding to a slight increase of the solid fraction C from 
0.785 to 0.79. Obviously, the small number of cycles ex- 
perienced in the present work does not allow to observe 
a saturation of the volume and eventually a critical state 
It is known that in similar cyclic solicitations [s^ . 
as well as in other experimental set-ups [lol 0| , even af- 
ter a large number of cycles (10"* cycles), the volume still 
slowly decreases. 

For slopes 9 < 6r, the granular pile is contracting 
: ey decreases. On the contrary, for slopes in the 
range [9r- 9a], the rearrangements cause the granular pile 
to dilate : ey increases. As a result, for small cy- 
cles {9rev = 10°) the compaction of the granular pile 
remains monotonous (Fig. IS)), whereas for large cycles 
i9rev = 18°), the granular pile exhibits compaction and 
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dilation stages (Fig. (S)). This latter behaviour is typical 
when shearing dense granular media, because shear de- 
formation is possible only if grains unjamm, confirming 
the onset of a coherent shearing above 9r. 

Interestingly, the overall behaviour of the pile over the 
complete solicitation is modified in the case of large cycles 
: the instantaneous compaction is twice more efficient 
for large cycles {dey/dt ~ —7 lO^^s"^) than for small 
ones (dcv/dt ~ —2 10^''^s^^). As a result and despite 
the dilation stage occurring when G [9j.;0a], the overall 
densification is larger for large cycles than for small ones 
for the same number of cycles, showing that densification 
depends on the solicitation path This leads us to 

speculate that grains are much more allowed to rearrange 
when the pile has previously evolved in the coherent shear 
regime identified prior to the unjamming transition. 
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FIG. 9: Kinetic energy Ek averaged over all the grains and 
slide-averaged over l°-intervals, normalised by the typical po- 
tential energy of a grain of mass m for a height d, as a function 
of 9 for 9rsv = 10° and 9rev = 18° (same symbols as figure0. 



This largest susceptibility for grains to rearrange in the 
case of large cycles {0rev — 18°) is confirmed by moni- 
toring the averaged kinetic energy of the grains Ek , rep- 
resented as a function of the pile slope in Figure |H| after 
slide-average over l°-intervals. The kinetic energy Ek 
involved in the local rearrangements, also occurring at 
small inclinations 9, is much higher when the granular 
pile has been loaded up above 0,- during its history. This 
evolution underlines the two different states and evolu- 
tions of the pile for the two cyclic amplitudes. 

Exploring the domain [9r',0a\, allows for stronger re- 
arrangements and larger volumetric strains during the 
full solicitation despite no apparent signature on ctt/cat. 
As a result, the volumetric deformation of the granular 
pile shows a strong dependence on the amplitude of the 
cycles. 



C. Evolution of the contact network 

We now investigate the effect of the deformations on 
the evolution of the micro-structure, i.e., the geometry 
of the contact network , th at are candidates at relevant 
internal variables i, H El IHlii El Isi . 
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FIG. 10: Coordinancy of the grains Z slide- averaged over 1°- 
intervals as a function of the slope of the pile for successive 
small and large cycles (same symbols as figure [7|l. 

The coordinancy Z, i.e., the mean number of contacts 
per grain, slide-averaged over l°-intervals is plotted in 
Figure [TUl as a function of the pile slope during the cyclic 
solicitations. For both values of the amplitude 9rev^ the 
coordinancy Z varies less than 1%. Nevertheless, it ex- 
hibits a weak hysteresis according to the sense of rotation 
: Z tends to decrease during loading, but to increase dur- 
ing unloading. 

Small and large cycles are identified by two different 
states : the value of the coordinancy depends on 9rev 
During small cycles {9 rev — 10°), Z remains approxi- 
mately at its initial value, by contrast with large cycles 
{9rev — 18°), for which a variation of the coordinancy 
occurs during the first exploration of the domain [9^.', 9a] 
and remains afterwards. 

Smallest coordinancy observed for 9rev = 18° corre- 
sponds to densest granular piles, contrary to intuitive 
observations [sij . Such an anti-correlation between the 
evolutions of the coordinancy and the volume under an 
external driving field could be related to the dependence 
of the packing equilibrium on its number of contacts. An 
"isostatic" packing has the minimal number of contacts 
required to remain at metastable equilibrium, by contrast 
with a "hyperstatic" one, that could loose some contacts 
and not its equilibrium. Yet for an initially static grain 
to rearrange, its contacts have to cooperate, that is easier 
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for a smaller number of contacts. In the hyperstatic pack- 
ing, a significant cooperation of contacts is needed for a 
grain to rearrange, by contrast with the isostatic pack- 
ing, in which a slight change in the force transmission at 
contacts will enable a local rearrangement of grains. In 
the case of large cycles, the significant decrease of Z ob- 
served during the first loading makes the packing become 
less hyperstatic : grains rearrangements occur easily, en- 
hancing the deformation of the pile. 



2.. Fabric of the pile 
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FIG. 11: Anisotropic intensity A (a) and fabric anisotropic 
direction $ of the contact network (b) as a function of the 
pile slope during the successive cyclic rotations for drm ~ 10° 
and 9rev = 18° (same symbols as figure [TJ. 

To further analyse the evolution of the micro-structure 
of the pile, we now consider the statistics of the orienta- 
tion of contacts. To do so, the fabric tensor t is computed 



using the following definition as defined in '525 : 

Q = l 

where i and j denote the space dimensions and a the con- 
tact, n is the unity vector normal to the contact surface, 
and n'^ is the total number of contacts. The fabric allows 
to investigate the evolution of the geometry of the con- 
tact network, especially if contacts are created, opened or 
modified such that their alignment tends towards a priv- 
ileged direction. In particular, the anisotropic intensity 
A is calculated as : 

A = 2x(ii-t2), (5) 

where ti and t2 are respectively the larger and smaller 
eigenvalues of t. In the case of an isotropic contact net- 
work, e.g., spatially periodical, the anisotropic intensity 
A = 0. The direction of the eigenvector of t correspond- 
ing to the largest eigenvalue is the fabric anisotropic di- 
rection, whose angle with the direction rCs is 

Figure lllh displays the evolution of the anisotropic 
intensity A as a function of the pile slope during the 
cyclic solicitations for 9rev — 10° and 9rev — 18°. Dur- 
ing small cycles, a rather isotropic state is observed, and 
only small variations of A with 9 are observed (A ~ 0.05) 
: the geometry of the contact network is only weakly af- 
fected by the solicitation and remains close to the initial 
one. On the contrary, for large cycles, the anisotropic 
intensity evolves significantly showing variations up to 
the maximum value A ~ 0.12. Approaching the unjam- 
ming transition obviously strongly modifies the geometry 
of the contact network, e.g., by introducing structural 
anisotropy. 

Figure lllt > displays the evolution of the fabric 
anisotropic direction $ as a function of the pile slope dur- 
ing the cyclic solicitations for 9rev = 10° and 9rev = 18°. 
Irrespective of the value of 9rev , evolves during the so- 
licitations and rotates in the sense opposite to the driving 
rotation : <i> decreases (resp. increases) when 9 increases 
(resp. decreases). Following this evolution, the fabric 
anisotropic direction tends to approach that of the free 
surface and for large inclinations, <I> tends towards the 
principal stress direction 4" = ±45°. This behaviour cor- 
responds to the mechanisms of creation of contacts in the 
direction of compressive stress, and the loss of contacts 
in the direction of extension. 

Besides these general observations, the contact net- 
work evolution exhibits more complex features. Both 
for small and large cycles, the fabric exhibits a strongly 
marked hysteresis, probing the irreversible modifications 
of the disordered micro-structure. The hysteresis loops 
of A and $ are very smooth. Furthermore when the 
granular pile has approached the unjamming transition, 
namely for large cycles {9rev = 18°), the shape of the 
hysteresis is dramatically changed, and the loop ampli- 
fied. This points out the peculiar effect of coherent shear- 
ing occurring prior to the unjamming transition on the 
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granular packing evolution. More particularly, the fabric 
evolves less rapidly during unloading than during loading 
for large cycles {9rev = 18°). In other words, the packing 
seems to be less able to reorganize its structure during 
unloading than during loading when having approached 
the unjamming transition in its history. The same re- 
mark applies both to the evolution of the anisotropic in- 
tensity A and to the fabric anisotropic direction $. This 
may be rather unexpected given that the efficiency of re- 
arrangements are increased in the case of large cycles, 
and this during the complete solicitation. 

Note that for large amplitudes {Orev = 18°), the hys- 
teretic cycles of A and $ are symmetric, by contrast with 
small amplitudes {9rev = 10°), for which an asymmetry 
related to the initial fabric anisotropy remains. Explor- 
ing the range of slopes [9r', Oa] when approaching the un- 
jamming transition allows a profound modification of the 
contact network according to its initial state, by increas- 
ing structural disorder. 

Altogether the fabric of the pile (see Z, A and $) at 
an angle 6, strongly depends on the solicitation path, as 
observed for the pile deformation (ey and Ek). 



V. RESPECTIVE ROLE OF STRONG AND 
WEAK CONTACTS 

Dense quasi-static granular media are known to exhibit 
two complementary contact networks, strong and weak, 
depending on the intensity of the normal force transmit- 
ted at contacts [H El El Hi ■ Very recently, the 
role of the critical contacts in the destabilization of a 
pile pri Is^ , as well as in the elasto-plastic response of 
a granular packing has also been underlined. Still, 
it was shown that during the slow tilting under gravity 
of a granular pile, critical contacts and more gene rally 
high frictional contacts are mainly weak contacts I2j.l33|. 
Here, in the fine of these studies [H 13 El El 
as a first step towards the understanding of the role of 
the micro-structure and the heterogeneous nature of con- 
tacts in the observed hysteresis, we choose to concentrate 
on the respective contributions of the strong and weak 
contacts to the hysteretic behaviour of the pile. Further 
works should definitely examine in more details the role 
of highly frictional contacts in the spirit of [H . 

Accordingly, the computation of the stress and fabric 
tensors are now restricted to strong or to weak contacts 
using respectively equation ^ or ^ In the present work, 
since a vertical gradient of the contact forces is related 
to the gravity field, a contact at a depth y is defined as 
strong (resp. weak) if it transmits a normal force larger 
(resp. smaller) than the averaged normal forces at depth 

y- 
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FIG. 12: Principal stress direction in the strong contact net- 
work as a function of the pile slope during the successive 
loading-unloading cycles for small and large amplitudes (same 
symbols as figure EJ. 
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FIG. 13; Fabric anisotropic direction of the strong contact 
network "1>^ as a function of the pile slope during successive 
loading-unloading cycles for small and large amplitudes (same 
symbols as figure 0. 



A. Strong contact network 

Figure El displays the principal stress direction in the 
strong contact network '^'^ as a function of the surface 
slope during small {6rev = 10°) and large {9rev = 18°) 
cycles. The evolution of is weakly hysteretic as a func- 
tion of the pile slope and the hysteresis is approximately 
the same irrespective of the amplitude of the solicita- 
tion, as observed previously for the total contact network 
(Fig. [7()) . The strong contact network represents by far 
the largest contribution to the global stress 0| . Accord- 
ingly and ^' exhibit very similar behaviours. 

The fabric anisotropic direction of the strong contact 
network is displayed as a function of the pile slope 
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during the cyclic solicitations in Figure [T^ The evolu- 
tion of exhibits hysteresis during the cycles of both 
amplitudes. However, unlike the previous observations 
on the fabric of the total contact network (Fig. \llb ) , the 
size of the hysteresis loops is small for drev = 10° and 
18°. The evolution of is very much correlated with 
the evolution of 'i'^, and the shape of the hysteresis loops 
of both and <f>* are smooth. 

During cyclic solicitations, neither the stress nor the 
fabric of the strong contact network is sensitive to the 
amplitude of the cycles. 



B. Weak contact network 



s 




FIG. 14: Principal stress direction in the weak contact net- 
work as a function of the pile slope during the cycles for 
small and large amplitudes (same symbols as figure EJ. 
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FIG. 15: Fabric anisotropic direction of the weak contact net- 
work as a function of the pile slope 8 during small and 
large cycles (same symbols as figure 0. 



Figure ^] displays the principal stress direction in the 
weak contact network vj/™ as a function of the surface 
slope during small {9rev = 10°) and large {9rev — 18°) 
cycles. The evolution of ^f"" is strongly influenced by the 
amplitude of the solicitations. This qualitative change 
is not reflected by the total contact network response 
(Fig. [7J)) due to the dominant contribution of the strong 
contact network to the global stress 0| . 

Figure El displays the evolution of the fabric 
anisotropic direction of the weak contact network <i>™ as a 
function of the pile slope for small {6rev = 10°) and large 
{drev = 18°) cycles. The response of the weak contact 
network is again strongly dependent on the cycles am- 
plitude 9rev , by contrast with the response of the strong 
contact network (Fig. I13|l . For small cycles {Orev — 10°), 
$™ exhibits nearly no hysteresis, and remains approxi- 
mately always normal to as observed in the case of 
continuous loading ^3§\ . On the contrary, for large cycles 
{9rev = 18°), a strong hysteretic behaviour is observed, 
which consists of a premature rotation of upon rever- 
sal : $™ switches beforehand during unloading in com- 
parison with loading. This effect leads to a remarkable 
orientation of the fabric of the two contact networks : $™ 
and are equal during a part of the cycle, attesting the 
temporary alignment of the privileged direction of the 
two contact networks. The particularly early evolution 
of the weak contact network geometry when the pile has 
approached the unjamming transition, is actually respon- 
sible for the delay observed in the reorganization of the 
total contact network. As a consequence, upon reversal, 
the packing keeps the memory of the previous orienta- 
tion of the total contact network during a large part of 
the cycle (Fig. El')- 

Contrary to the smooth shape of the hysteresis loops in 
the strong contact network (Fig. I13|l and in the total con- 
tact network (Fig. II lb ') , the hysteresis loops in the weak 
contact network (Fig. I15|) exhibit an abrupt macroscopic 
jump. It should also be mentioned that the weak contact 
network response exhibits large fluctuations from one re- 
alization to another one. These fluctuations are strongly 
enhanced close to the unjamming transition. Evidences 
of these fluctuations are for instance the noisy signal for 
(Pig |i4|) and the absence of superposition of the hys- 
teretic loops for (f)"" (Fig. El) ■ 

Altogether, in the weak contact network by contrast 
with the strong one, both stress and fabric are very sen- 
sitive to the amplitude of the cyclic solicitations. Two 
conclusions can be drawn from these observations : 

• In the limit of small solicitations, e.g., for cycles 
with 9rev = 10°, the hysteretic behaviour of the 
granular pile when looking at both stress and fabric 
state is dominated by the strong contact network 
contribution ; 

• For large cycles, with Orev > Or, namely when 
the pile approaches the unjamming transition, the 
response of the weak contact network becomes 
strongly affected, and the hysteretic behaviour of 
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the pile is strongly modified despite the lack of 
modification of the response of the strong contact 
network. 

Accordingly the hysteretic behaviour of the pile is the re- 
sult of a complex interplay between the strong and weak 
contact networks, whose relative contributions strongly 
depend on the amplitude of the cycles. 

VI. DISCUSSION AND CONCLUSION 

The successive tiltings of the granular pile from dif- 
ferent initial conditions up towards the avalanche angle 
allow to demonstrate that the density and the correla- 
tion length of the critical contacts are relevant dynamic 
response functions to describe the metastable evolution 
of the packing towards the unjamming transition under 
external loading. Due to the transient nature of the crit- 
ical contact state, these quantities are representative of 
actual loading of the pile, but they do not characterize 
its plastic state. 

The response of the granular packing to the quasi- 
static cyclic solicitations reveals that the evolution of 
the critical contacts is hysteretic when the pile has ap- 
proached the unjamming transition. At the transition 
from loading to unloading, critical contacts can not any 
more sustain the sliding condition, so that the frictional 
state of the packing is immediately removed, attesting 
the end of coherent shearin g, a s observed in the case of 
experimental shear reversal |2.1L l37| . 

When looking for such memory effects on other fea- 
tures of the pile, namely on stress, strain and fabric, they 
exhibit different responses. The stress state exhibits weak 
hysteretic evolutions as a function of the slope, and the 
hysteresis reveals to be independent of the amplitude of 
the cycles. On the contrary, the strain and the fabric evo- 
lutions of the pile, related to structural disorder, are very 
sensitive to the amplitude of the cycles. These observa- 
tions will hardly be explained by a simple constitutive 
law relating only stress and strain of the granular pile : 
such a relation should be not so trivial, but should take 
into account the anisotropy of the pile micro-structure 
EE El E E Hi Ei . The observation of piles with the 
same stress state but characterized by different granular 
fabrics, would suggest that stress induced by gravity are 
transmitted in the pile whatever the geometry of the to- 
tal contact network, and whatever the structural disorder 
level. 

Further analyses of the strong and weak contacts con- 
tributions to the macroscopic behaviour allow to clarify 
these intriguing observations. In the strong contact net- 
work, neither the stress nor the fabric depends on the so- 
licitation amplitude. In the weak contact network, both 
the stress and the fabric are sensitive to the cyclic ampli- 
tude. The apparent mismatch between stress and fabric 
at the scale of the total contact network does not sur- 
vive when separating weak and strong contact networks 
contributions. It simply comes from the strong contact 



network weight in the overall stress [3^ . Accordingly, the 
above analyses stress the relevance of a constitutive law 
which specifies the two-phases nature -weak and strong- 
of granular media. In the strong contact network, stress 
and fabric are very much synchronised, and a rather sim- 
ple relation could relate them. The behaviour of the 
weak contact network is much more complex. A strong 
memory of the unjamming transition affects the response 
of the weak contact network during the full solicitation, 
probably due to high sensibility of weak contacts to co- 
herent shearing occurring prior to the unjamming transi- 
tion. Finally the observed fluctuations in the weak con- 
tact network are presumably related to long-range cor- 
relations. These correlations within the pile allow a lo- 
cal rearrangement of grains to further propagate and to 
affect essentially the weak contacts, due to their high- 
est frictional mobilization, as recently shown in l33| . 
This high sensibility of the weak contacts to local pertur- 
bations enables them to evolve collectively and to change 
the micro-structure of the packing. A signature of the lo- 
calization of the correlations in the weak contact network 
could be seen in the particular shape of the hysteresis 
loops. Whereas the loops are very smooth in the total 
and strong contact networks, the hysteresis loops in the 
weak contact network exhibit abrupt macroscopic jumps. 
These respective delay and instantaneous responses could 
be related to the relative contribution of structural dis- 
order and finite-size interactions or correlations, as sug- 
gested in [s^ |53|. In this picture, the hysteresis with 
a jump in the weak contact network would suggest that 
the size of correlations is larger than the typical length- 
scale of the structural disorder in this phase, contrary to 
the predominant role of disorder in the micro-structure 
of the packing. Altogether, such complex behaviours - 
spatial correlations and memory effects- call for more 
elaborated constitutive laws as in [i lliTllisL ligllsollssll . 

To conclude, the behaviour of the pile is strongly mod- 
ified close to the unjamming transition, as experimen- 
tally observed |23 . l2a . l3l] | , coinciding with important ir- 
reversible modifications of the micro-structure. These 
modifications result from a specific solicitation of weak 
contacts : these latter play a considerable role in the be- 
haviour of the pile despite their marginal contribution 
to the stress state. The present paper demonstrates the 
relevance of a two-phases description, not only for the 
destabilization process, but also for the quasi-static rhe- 
ology, as soon as the granular sample has approached the 
unjamming transition in its past history. Further stud- 
ies should investigate the complex behaviour of the weak 
phase in order to extract a constitutive relation for gran- 
ular materials, accounting for their two-phases nature. 
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